Discrete solitons and vortices in hexagonal and honeycomb lattices: 
existence, stability and dynamics 
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We consider a prototypical dynamical lattice model, namely the discrete nonlinear Schrodinger 
equation on non-square lattice geometries. We present a systematic classification of the solutions 
that arise in principal six-lattice-site and three-lattice-site contours in the form of both discrete multi- 
pole solitons and discrete vortices. Additionally to identifying the possible states, we analytically 
track their linear stability both qualitatively and quantitatively. We find that among the six-site 
' configurations, the "hexapole" of alternating phases (0-7r), as well as the vortex of topological charge 

5* = 2 have intervals of stability; among three-site states, only the vortex of topological charge 
, S = 1 may be stable in the case of focusing nonlinearity. These conclusions are confirmed both for 

hexagonal and for honeycomb lattices by means of detailed numerical bifurcation analysis of the 
stationary states from the anti-continuum limit, and by direct simulations to monitor the dynamical 
instabilities, when the latter arise. The dynamics reveal a wealth of nonlinear behavior resulting 
not only in single site solitary waveforms, but also in robust multi-site breathing structures. 



I. INTRODUCTION 



Hamiltonian lattice or quasi-discrete systems have become popular in the last few years, to a considerable extent 
due to experimental implementations of such systems drawn from various branches of physics. One of the first 
examples where such developments became relevant was in the nonlinear optics of fabricated AlGaAs waveguide 
arrays There, the interplay of inherent discreteness and nonlinearity led to the emergence of numerous interesting 
phenomena including Peierls-Nabarro potential barriers, diffraction and diffraction management @, gap solitons 0], 
and so on (see also the reviews [1, Q and references therein) . 

On the other hand, more recently another area of nonlinear optics that has been central to the development of 
both theoretical as well as of computational tools to study such systems, has been the setting of optically induced 
^vq , photonic lattices in photorefractive crystals such as SBN. There, the original theoretical proposal of these lattices 
Q was susbequently followed by experimental realizations @, \f^, paving the way for the observation of a diverse 
array of novel and interesting phenomena in such crystals. These include the formation of patterns such as dipole 
@, quadrupole [1(3] and necklace [ll[ solitons, impurity modes discrete vortices [HI, [lij], rotary solitons flEl ]. 
higher order Bloch modes |lo| and gap vortices [13], the observation of two-dimensional (2D) Bloch oscillations and 
Landau-Zener tunneling Il8l . the observation of localization and diffraction in honeycomb [lj5[ , hexagonal [2(| and 
quasi-crystalline lattices |2l| . and most recently the stud y o f Anderson localization in disordered photonic lattices [22[ 
(for a review of some of this activity see, e.g., Refs. [23|,[24j]). 

Finally, we should note that similar dynamical lattices have also become of interest in an entirely different area of 
physics, namely the atomic physics of Bose-Einstein condensates (BECs), when trapped in periodic potentials; see, 
e.g., the recent reviews [25l. [26l. |27| . 

In this work, we will focus, in particular, on the recently e merg ing area of non-square lattices in waveguide arrays, 
as well as in light induced photonic crystals, [H, US IH> HI; H3]- Such lattices were also considered earlier from a 
theoretical perspective in discrete settings corresponding to waveguide arrays [3l|. We will study here the existence, 
stability and dynamical properties of multi-pulse solitary wave structures, as well as of discrete vortex structures in 
both hexagonal and honeycomb lattices. We will focus on two prototypical contours of such lattices. Namely, a more 
extended six-site contour, as well as a reduced three-site contour, both depicted in Fig. [T] Our prototypical model 
of interest will be the discrete nonlinear Schrodinger (DNLS) equation and our results will be presented for the case 
of a focusing nonlinearity; however, our findings can be directly transformed to the case of a defocusing nonlinearity. 
Additionally, we should note that similar results can be obtained in Klein-Gordon models and have been illustrated, 
e.g., for three-site contours in hexagonal lattices [32| . 

It is relevant to note that crystalline configurations of strongly coupled doped plasmas (dusty_plasma crystals) 
occur in the form of ID or 2D monolayers formed in low-temperature gas discharge experiments [331 ]. Interestingly, 
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such dust crystals generically appear as spontaneously formed hexagonal 2D arrangements [34[ , although alternative 
configurations also include honeycomb 2D lattices [35| and ID dust chains (when appropriate trapping potentials 
are used for lateral confinement j36|). A discrete Klein-Gordon description has recentl y b een employed to model the 
dynamics of transverse vibrations of dust grains in dusty plasma crystals, both in ID [37| and in hexagonal 2D dust 
lattices [38j]. 

The key findings that we report here are the following: 

• For the focusing nonlinearities considered herein, in a six-site honeycomb /hexagonal contour, topological charge 
S = 2 configurations may be stable, while S = 1 ones can never be stable. This represents a notable qualitative 
difference from the results in the case of a square lattice [39| , where the prototypical contour consisting of four 
sites features a potentially stable S = 1 vortex, as well as a genuine potentially stable analog of the S — 2 vortex 
(the so-called S — 2 quasi- vortex (40l|). 

• In these contours, also six-site alternating 0-7r phase configurations are potentially stable, while in-phase config- 
urations are not stable. While this instability can be implicitly inferred from the instability of the corresponding 
building blocks (i.e., the instability of the in-phase dipole and the potential stability of the out-of-phase dipole 
41]), its quantitative characteristics can only be traced through the approach presented below. 

• In three-site contours, the only potentially stable configuration is that of a discrete vortex, while both in-phase 
and alternating phase configurations are observed to be unstable. 

• The evolution of the dynamical instability in these lattices is more complex than in the square lattice case, and 
may involve not only degeneration to single-site solitons but possibly to multi-site solitary wave structures, and 
often the formation of robust breathing states, consisting of multiple sites (possibly even as many as the original 
configuration). In fact two clear breather formations recur in multiple simulations: 

— Two sites with fluctuating, usually opposite, relative phases and oscillating amplitudes of comparable 
magnitude. 

— Two sites with different amplitudes oscillating between the same relative phases and opposite relative 
phases, depending on whether the amplitudes are closer or further, respectively. 

Six-site configurations with phases of or it will be collectively called "hexapoles" herein, while three-site configura- 
tions with phases or tt will be collectively termed "tripoles". 

Our presentation of the above findings is structured as follows. In section II, we provide the background for 
the theoretical analysis. In section III, we corroborate the theoretical findings with numerical bifurcation results 
illustrating the various nonlinear modes in both hexagonal and honeycomb geometries and their stability properties, 
while the corresponding dynamics are presented in section IV. Finally, in section V, we summarize our findings and 
present our conclusions. 

II. THEORETICAL ANALYSIS 

We consider the following discrete nonlinear Schrodinger equation (DNLS) for the two geometries of interest: 

sA d - I (i) 



.dllrn,n A 1 



dz 

U m ',n> - \N\u m ,n | , (2) 

where the summation is over the set iV of nearest neighbors (denoted by (m',n')) of the site (m,n), \N\ is the 
cardinality of the set N of neighboring sites (six in the hexagonal geometry and three in the honeycomb), and u m n 
models, e.g., the envelope of the electric field in the corresponding waveguide [3l|, while e represents the coupling 
strength between nearest neighbor nodes; z denotes the propagation distance along the crystal. 

In the, so-called, anti-continuum (AC) limit e — > the sites are uncoupled. In this case, explicit solutions over 
contours of nodes indexed by j can be easily found in the general form Uj = \/~Kexp(i6j) exp(iAz), where A is the 
propagation constant and 8j 6 [0, 2ir). Without loss of generality we will fix the value of A, taking A = 1, and we will 
consider three- and six-site contours in each of the hexagonal and honeycomb geometries, shown in Fig. [TJ According 



FIG. 1: (Color online) Discrete lattice configurations for the hexagonal geometry (left), in which each node has six neighbors, 
and the honeycomb geometry (right), in which each node has three neighbors. The relevant "hexapole" configurations are 
represented by the red circles and the "tripoles" are given by blue squares. Notice that the relevant three site configuration for 
the honeycomb is composed of next-nearest neighboring sites. 



to the earlier work of [39|, E2J (see also [43J), the necessary condition for a discrete contour M in the AC limit to 
persist in the presence of non-zero coupling is: 

g, = sin(0j - 9 j+ x) + sin(0,- - = 0, (3) 

for all j € M. The stability can also be determined from the eigenvalues jj of the \M\ x \M\ Jacobian Jjk = dgj/d9k- 
In particular, for each eigenvalue jj, the full linearization of Eq. |2]) around a stationary solution with non-zero nodes 
in M will have eigenvalue pairs Xj given, to leading order, by the following relation in the case that the sites in M 
are nearest neighbors: 



X 3 =±^~e. (4) 

If the non-zero sites comprising the contour are next-nearest neighbors instead, as in the case of the three site 
contours for the honeycomb lattice geometry (see Fig. [T]), then e is replaced by e 2 in the previous relation. Unstable 
solutions for weak coupling (small e) can then be identified as those for which the eigenvalues Xj have non-zero real 
part, given the Hamiltonian nature of the model. The Jacobian matrix has the following form: 

7j) + cos(9j-i-9j), j = k, 
(J)j,k= { -cos^-flfe), j = k±l, (5) 

\k-j\>2. 

We will consider primarily contours M such that \6j + i — 9j\ — A6 is constant for all j £ M, \9i — = A9 and 
A9\M | = mod 2ir, except one case which will be treated separately. In the primary case, all the non-zero elements 
of this matrix are then factors of a — cos(A#), and the eigenvalue problem of the Jacobian J reduces to the following 
difference equations: 

a{2x n - x n+ i - z n _i) = 7jX„. (6) 

These can be solved by a discrete Fourier transform with any eigenvector x n ~ exp(i2imj /\M\), whence jj — 
4asin 2 (7rj7|M|) and then 




= ±y 8ecos (A9) sin 2 ^j^jj ■ (7) 

Recall that for the honeycomb three-site next-nearest-neighbor contours the above formula should be used with 
e replaced by e 2 . The special case of the three-node contour with phases 0, it and can be treated also in the 
framework of the Jacobian of Eq. ([5]) [and its eigenvalues computed by Eq. Q], although it does not fall under 
the general calculation of Eqs. ©-([7]). We will consider nodes separated by either A9 = 0,7r,7r/3, or 27r/3, for the 
different contours in this work. 

These predictions for the linear stability eigenvalues will be compared to numerical results for the linear stability 
of the stationary solution i> mi „ exp(iz) of Eq. The stability will be analyzed by using the following ansatz, 

u m ,n = e iz [v min + S (p m , n e Xz + 9™,„e A * z )] , (8) 
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FIG. 2: (Color online) Six-site real-valued configurations in a hexagonal geometry. The top row corresponds to the unstable 
A9 — 0, or "in phase" solutions, while the bottom corresponds to the stable A9 = n, or "out of phase" ones. From left, the first 
column is the profile at e = 0.08, the second column shows the corresponding linearization spectrum (X r ,Xi) of the eigenvalues 
A = A r + iXi, and finally the third column shows the continuation in e of the actual eigenvalues [real, A r , and imaginary, A;, 
components] (solid) and the theoretical predictions given by Eq. (dashed). 



and solving the ensuing eigenvalue problem for the eigenvalue A and the eigenvector [p m ,n, Qm.n] T ■ In the above, the 
asterisk (*) denotes complex conjugate and "T" denotes transpose. Notice that when the resulting eigenvalues of the 
linearization possess a nonzero real part, the solution will be exponentially unstable, with a growth rate corresponding 
to the real part of the most unstable eigenvalue. 



III. NUMERICAL RESULTS 



Both in this section, detailing the various configurations and their corresponding stability over the six-node and 
three-node contours, and in the next one, comparing the corresponding dynamics, we will partition our discussion to 
two subsections. The first one will be devoted to the results obtained for the hexagonal geometry, and the second 
devoted to the case of the honeycomb geometry. 



A. Hexagonal geometry 

First, we will study six-site contours, of which we will consider four. The first two of these are real and are such 
that either AO = or A9 = n (any additional combination of and 7r phases is also possible but the main qualitative 
characteristics of stability will not change from those reported below). The relation ([7]) for AO = predicts double 
eigenvalue pairs at ±\/2e and ±vfe and single pairs at ±y/8e and 0, while for AO = it each of these is multiplied 
by the imaginary unity. Direct numerical computation and continuation in the coupling parameter e from the AC 
limit confirm the predictions presented in Fig. [2] The in-phase configuration with AO = becomes strongly unstable 
away from the AC limit, while the out-of-phase configuration with AO — it is stable for small e. It should be noted 
that more generally any configuration that has two adjacent in-phase nodes along the six-site contour will also be 
unstable for all values of e, while the only potentially stable configuration of this type (real solution comprising and 
7r phases) is the out-of-phase adjacent node structure of AO — tt. However, even for that configuration, the imaginary 
eigenvalues which bifurcate from the origin in the AC limit have the topological property of, so-called, negative Krein 
signature [39]; this means practically that they become structurally unstable upon collision with other eigenvalues, 
such as those of the phonon band, which have positive Krein signature. Hence, when the coupling becomes sufficiently 
large (e > 0.06), these eigenvalues eventually intersect with the continuous spectrum (the phonon band) edge located 
at ±zA, and result in Hamiltonian-Hopf bifurcations associated with complex quartets of eigenvalues and oscillatory 
instabilities. Such collisions can be detected in the graphs illustrating the lowest imaginary eigenvalue parts Xi in the 
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FIG. 3: (Color online) The six-site vortices in the hexagonal geometry are shown. The left four images show the unstable singly- 
charged solution branch with A9 = 7r/3, while the right set is for the stable doubly-charged solution branch with AO = 2n/3. 
The singly-charged vortex is unstable, while the doubly-charged one is stable (until the oscillatory instability resulting from 
the collision of the pairs of negative Krein signature with the phonon band). The top row of each set displays the modulus 
(left) and argument (right) of the solution with e = 0.025 (top) and e = 0.125 (bottom), while the bottom left is the linearized 
spectral plane and the bottom right is the continuation of the relevant eigenvalues from the AC limit, with the solid and dashed 
lines representing the numerical solution and the theoretical prediction, respectively. 



bottom right panel of Figured and are associated with the points where the eigenvalue trajectories begin to level out. 
The spectrum of the solution at e = 0.08 in the bottom middle panel of Fig. [5] reveals the presence of such quartets. 

Next, we consider the complex valued solutions along the six site contours, for which our conditions guarantee 
vorticity (i.e., the relevant solutions will be discrete vortices whose phase completes a round trip of a multiple of 2tt 
along the discrete contour). The fundamental solutions here are for AO — 7r/3, which is a singly-charged vortex, and 
A8 = 2tt/3, which is a doubly-charged vortex. The relation |(7J) predicts that the singly-charged vortex will be unstable 
with double eigenvalue pairs ±y/e and ±y/3e, and single pairs at ±V4s and 0. On the other hand, the doubly-charged 
vortex is actually stable in this lattice geometry for sufficiently small values of the coupling, with the same pairs as 
the singly charged vortex, except multiplied by i. Figure [3] presents both types of configurations, indeed illustrating 
the numerical linear instability of the former, and numerical linear stability of the latter structures. Nevertheless, 
it should be pointed out that for higher values of the inter-site coupling (e > 0.1) in this case also, the topological 
charge S — 2 solution eventually becomes unstable as well due to oscillatory instabilities, as is shown in the bottom 
right continuations of the relevant eigenvalues of Figure [3J Notice also the generally excellent qualitative and good 
quantitative agreement -at least for small values of e (for larger values of the coupling parameter higher order effects 
become important)- between the theoretical predictions of Eq. ([7]) and the numerical results. 

We now turn to the configurations comprised of three lattice sites. We consider three such cases, similarly to [32| 
where a Klein-Gordon model was considered. The first two are the standard real (A6 = 0) and complex-valued 
(A8 = 27r/3, corresponding to a discrete vortex of topological charge S = 1) ones, while the last one is the non- 
standard case with phases 0, ir and for the three sites. For A0 — the theoretically predicted double pair of ±v6e 
and pair at are confirmed to exist in the middle row of Figure [H while for the discrete singly-charged vortex solution 
with A9 = 2tt/3, the double pair ±\/3si and a pair at are also found to reasonably approximate its linearization 
eigenvalues in the bottom row of Fig. [4] It should be noted that for larger coupling (e > 0.02) the double pair of 
eigenvalues splits, as can be observed in the numerical results. These solutions also become unstable as usual after 
the first eigenvalue collides with the continuous spectrum at e » 0.09. On the other hand, for the configuration with 
phases 0, n and of the top rows, one pair of eigenvalues remains at the origin, due to the phase invariance, while 
one of the remaining two pairs becomes imaginary (theoretically predicted as iiVfe) and the other one becomes real 
(predicted as ±v2~i)- We can again observe that the full numerical results agree well not only qualitatively but also 
even quantitatively with the theoretical description. 
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B. Honeycomb Geometry 

We now explore the same configurations systematically in the case of the honeycomb lattice geometry, in which 
each node has three neighbors as opposed to six. In. Fig. [5J we again consider two representative real configurations, 
namely the in-phase six-site structure (top row), and the out-of-phase hexapole, where adjacent neighbors have a 
relative phase shift of %. We find that the principal stability characteristics are similar to those in the hexagonal 
case, in agreement with the theoretical prediction. In fact, we can observe that quantitatively the agreement of the 
linearization eigenvalues is arguably even better between theoretical predictions and full numerical computations in 
this case. This is because the central site, mediating second-order inter-site interaction between the excited sites, is 
absent in the six-site configuration in the honeycomb lattice (contrary to the case for the hexagonal configuration). 
This feature reduces the role of higher order corrections to the theoretical predictions and hence renders the leading 
order predictions accurate for wider parametric ranges. This is a feature that we consistently observe throughout our 
honeycomb lattice results. 

The six-site discrete vortices are illustrated for this lattice in Fig. [6l Once again, the interesting feature of the 
immediate and generic (i.e., independently of the precise value of e) instability of the vortex of topological charge 
5 = 1 can be observed, while the vortex of topological charge S = 2 is stable for small values of the coupling and is 
destabilized for sufficiently large couplings (e > 0.13) by means of oscillatory instabilities. Note the larger value of 
the coupling necessary for the onset of an oscillatory instability here as compared to the hexagonal case, presumably 
a result of the higher order terms present in the latter due to presence of the center site. 

Finally, the interesting feature of the three-site configurations in the honeycomb case is that they now constitute 
next-nearest-neighbor configurations. As a result, the theoretical prediction that should be compared to the full 
numerical results is now oc e rather than oc y/s. This is clearly seen to be consonant with the full numerical findings 
of Fig. [7J not only for the strongly unstable (with a double real pair iVfe) configuration of the in-phase case, or for 
the linearly stable (for e < 0.43)) vortex case (with a double imaginary pair ±\/2ie), but also for the top-row, 0-7T-0 
case of one real (±\/2e) pair and one imaginary pair of eigenvalues. 

IV. DYNAMICAL EVOLUTION RESULTS 

We now examine the nonlinear dynamics of an unstable solution of each configuration upon integration of a slightly 
perturbed waveform u = u s (l + u r ), where u s is the complex valued vector field which is a stationary unstable 
solution to Eq. and u r is a random noise field (i.e., a field in which every entry is a random variable distributed 
uniformly in the interval between ±0.05max{ m n }[|u mj „(t = 0)| 2 ]). Since the coupling sensitively affects the dynamics 
(in particular, larger coupling facilitates communication between sites and hence propagation of the instability), we 
use a fixed coupling of e = 0.1 for all solutions except for those which are stable until larger values of e. In a few 
seemingly counterintuitive cases we examine the cases of larger perturbation and coupling, and find that these effects 
(more so the coupling) do indeed influence the dynamical evolution. We will see that several cases degenerate to 
similar two-site breathing structures with phase correlation which may be either out-of-phase or oscillating between 
in- and out-of-phase. 

A. Hexagonal Geometry 

First, we explore the evolution of characteristic unstable solutions from the families of configurations in a hexagonal 
geometry given in section [III A[ Within this class we begin with the six-site configurations. The evolution of the real 
valued solution with A8 = from the family in the top row of Fig. [2] is displayed in Fig. [U The rapid destruction of 
the original configuration confirms the linear stability analysis, which predicts strong instability from five pairs of real 
eigenvalues. However, for e = 0.1 (top set) and a 5% perturbation, after the destruction of the initial configuration, 
a robust three site oscillating breather state emerges (note the plot of individual site amplitudes as a function of 
propagation distance in the third row of Fig. [8]). Despite the apparent coherence of the amplitude oscillations, the 
relative phases of the sites appear to be uncorrelated and are not shown. A similar phenomenon is observed for a much 
larger initial perturbation of 25% of the initial amplitude (bottom left panels of Fig. [5]), although here the amplitude 
oscillations remain irregular even with three populated sites, and after a long distance, a nonlinear dynamical structure 
emerges in the form of a two-site breather. Again, however, there is no definite pattern in their relative phases. For 
a much larger coupling, on the other hand, as shown in the bottom right panels, all sites except for one decay very 
rapidly and a single site survives for long distances. Dynamical evolution of a real valued solution from the bottom row 
of Fig. [2] with A9 = ir is displayed in Fig. [9] The original configuration takes considerably longer to decompose than 
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the in-phase counterpart given above, confirming the expectation based on the small magnitude complex quartet of 
unstable eigenvalues of the linearized system. Once again, in this case a two-site structure with oscillating amplitudes 
persists long after the original break-up. However, in this case, there is a strong phase correlation, and when the 
amplitudes of these sites are close they are in-phase, while when they are distant they are out-of-phase (shown in the 
right panels). 

Next, we consider the vortex solutions with six sites. Both of these configurations confirm again the linear stability 
analysis, and also both feature two-site breathers for long distances. The singly-charged vortex (AO = tt/3) from the 
left panels of Fig. [3] decays into a breather with uncorrelated phases, similarly to the bottom left panel of Fig. [5] and, 
hence is not shown. The evolution of the more stable doubly-charged vortex (AO = 27r/3) from the right panels of 
Fig. [3]is in Fig. [TU] Notice the almost harmonic oscillations of the breather shown in the inset for the doubly-charged 
case. Here the two sites are also of comparable amplitudes, but as they oscillate they remain usually out-of-phase 
with each other as shown in the right panels. Another feature of both of these cases is that one of the two ultimately 
surviving sites is the originally unpopulated center site, which inherits mass from other sites when they decay (see 
also the insets in each figure). 

We now consider the three-site configurations. Both the Oi — 0, 7r, solution from the top rows of Fig. |4] (unstable 
due to one real pair of eigenvalues) and the more unstable AO = solution given below that (which is unstable due to 
two pairs of real eigenvalues), ultimately decay into in-phase/out-of-phase breathers such as the one shown in Fig. [9l 
although faster in the latter case due to the two unstable directions. The latter is displayed in Fig. QT] The final three 
site configuration is from the potentially stable AO = 2tt/3 family in the bottom panels of Fig. 0J The imaginary 
eigenvalues with negative Krein signature do not reach the continuous spectrum until a large coupling value in this 
case, and so we investigated the dynamics for e = 0.2. Despite the magnitude of the growth rate being comparable 
with the previous cases, two of the original populated sites here rapidly decay and a robust single site remains. This 
may be a result of the stronger site interaction induced by the larger coupling. 

B. Honeycomb Geometry 

We now turn to the same configurations as above but in the honeycomb geometry, as explored in section IIIIBI 
Interestingly, in this case, for e = 0.1, and a 5% perturbation, all configurations result in multi-site breathing 
structures with up to four populated sites for long propagation distances. Since the dynamical evolution of the 
six-site configurations in the hexagonal geometry all involve communications with the center site, it is reasonable to 
hypothesize that this is a major contributor to the rather significant differences observed below between the nonlinear 
evolution presented in this and the previous subsection. 

First, we display the results of the evolution of a real valued configuration from Fig. [5] in Fig. [T2J The linearized 
system of the solution with AO = in Fig. [TJ] is strongly unstable with five real pairs of eigenvalues, and the one 
with AO = 7r has all the same multiplied by the imaginary unity. The dynamical evolution confirms the stability 
analysis and two sites decay very rapidly for the in-phase configuration, while much more slowly for the more stable 
out-of-phase one (not shown). On the other hand, it is noteworthy that the four sites persist for long distances in 
each case and that the more linearly stable out-of-phase one decays into three sites eventually, which have apparently 
uncorrelated phases. The resulting four site breather in the in-phase case for e = 0.1 (left) is actually comprised 
of two out-of-phase breather pairs, such as the one in Fig. \TU\ while for e = 0.3 (right) the phases of the unequal 
amplitude breather pair oscillate between in-phase and out-of-phase. Also, as in the hexagonal case of the AO = 
family, we explored the sensitivity of the nonlinear evolution to larger perturbation and coupling and found that two 
sites robustly remain for the larger coupling e = 0.3, while four remain for e = 0.1 (not shown). For this reason, these 
solutions were continued for an extra long propagation distance up to z = 2000, and for consistency and comparison 
the remaining cases in this setting will also be continued for the same distances. 

The instability of the discrete vortices from Fig. [S]results in multiple sites persisting with large amplitude oscillations 
for long distances, ultimately evolving to an out-of-phase two-site breathing structure for the singly-charged (AO — 
tt/3) one (not shown) and a four-site structure for the doubly-charged (AO = 27r/3), shown in Fig. [T3J This four 
site structure consists of an out-of-phase pair close in amplitude and an unequal amplitude pair oscillating between 
in-phase and out-of-phase (see right panels) until a very long distance when they reshape into two out-of-phase pairs 
(phase not shown). The latter part is reminiscent of the result of evolution of the in-phase hexapolc for small e given 
in the left panels of Fig. [TJJ 

Finally, we show the evolutions of the three-site configurations from Fig. [7] Figure [T3] displays the dynamics of 
an unstable 0i = 0,7r, solution. The persistence for very long distances of the three sites for the smaller coupling 
prompted an investigation of a solution with larger coupling of e = 0.27 from this family. This turned out to decay 
very rapidly to a single site (not shown). In the smaller coupling case, an intricate breathing pattern emerges which 
apparently converts the mode into a three-site breather (rather than a three-site stationary solution) . The three-site 
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configuration with AO = is not shown, but again here all three sites survive for a long propagation distance for 
e = 0.1. This does not necessarily contradict the linear instability, since the configuration deviates almost immediately 
in terms of amplitude distribution. There is no clear correlation in the phases in this case. Again the persistence 
of all three sites for e — 0.1 prompted investigation for a larger coupling e — 0.3 and again a single site ultimately 
remained, although in this case two sites also persisted for a significant distance before the ultimate degeneration 
into a single-site waveform. For the last three site configuration the same consideration of the coupling arises, since 
this configuration is unique among those considered here, in the sense that a considerably larger coupling strength 
is required for the imaginary eigenvalues to collide with the continuous spectrum and the instability of this state to 
occur. Even with the very mild instability when the first imaginary pair collides with the phonon band at the very 
large coupling value of e = 0.43, the dynamics clearly illustrate the oscillatory instability. The original configuration 
persists until z — 30, ultimately concentrating primarly on a single site for long propagation distances. Aside from the 
six-site in-phase configuration, this is the only one for which the dynamics are qualitatively similar in the honeycomb 
and hexagonal geometries. Breathers remained for certain parameter values for all other configurations considered. 
The relative phases of the two-site breathers which recurred in many of the simulations suggest that these may exist 
as potentially stable time-periodic solutions. 

V. CONCLUSIONS 

In the present work we have examined both discrete soliton and discrete vortex configurations on non-square lattices 
for a prototypical discrete nonlinear Hamiltonian evolution equation (namely, the DNLS equation) of diverse interest 
to various areas of nonlinear optics and potentially of atomic physics. We studied, in particular, three- and six-site 
configurations in non-square lattice geometries in the case that each node has six neighbors (hexagonal lattice), as 
well as in the case that each node has three neighbors (honeycomb lattice) . Theoretical predictions for the stability of 
six-site real configurations were the same in each case and the in-phase version was strongly unstable (and generally all 
configurations with any pair of in-phase nearest neighbors will be unstable) , while the out-of-phase configuration was 
subject only to weak oscillatory instabilities stemming from complex quartets of eigenvalues that arise when imaginary 
eigenvalues with negative Krein signature collide with the continuous spectrum (of positive Krein signature). On the 
other hand, among complex solutions, we highlighted the cases of topological, singly- and doubly-charged structures 
(discrete vortices). In both the hexagonal and honeycomb geometries it was found that the former configuration is 
strongly unstable, while the latter may be stable for sufficiently weak couplings. It is also relevant in this context to 
point out that our results were presented for the case of a focusing nonlinearity, but it is straightforward to extend 
them to the defocusing case. There, it is expected that these features will be reversed, i.e., the vortex of charge 
5=1 will be stable, while that of S — 2 will be unstable in contrast to the effect of such a transformation on 
the fundamental four-site contour in a square lattice. Similarly the in-phase structure will be stabilized, while the 
out-of-phase one will be destabilized which, however, is consistent with the case of the square lattice. This can be 
inferred by a staggering transformation along the one-dimensional contour of excited sites, which changes by tt the 
phase of, say, just the odd (or equivalently just the even) sites of the contour. For the three-site configurations, on 
the other hand, there was a difference between the honeycomb and hexagonal geometries which arises due to the 
fact that such a configuration occurs in the honeycomb case comprised of next nearest neighbors and therefore the 
eigenvalues grow at a higher order in e: in fact, linearly (compared to a growth proportional to \fe in the hexagonal 
case). Nevertheless, in all cases, our theoretical predictions were confirmed qualitatively, but also quantitatively (at 
least for small e) by numerical continuation/bifurcation results. 

Finally, the dynamics of these states revealed a number of somewhat unexpected features. Unlike the expectation 
of square geometries that the relevant configurations may typically degenerate to a single-site solitary excitation, 
we have observed in various cases, that not only can the configuration reduce to waveforms with a larger number 
of excited sites, but also it can even preserve its original amplitude profile in terms of the number of excited sites 
- exhibiting breathing oscillations in the amplitude of these sites, or an internal reshaping of their relative phase. 
A breather composed of two-sites having comparable amplitudes and being usually out-of-phase with one another 
recurs in the dynamics of several solutions, as well as one with two uneven amplitudes in which the phases oscillate 
between in-phase and out-of-phase according to whether the amplitudes are closer or further apart, respectively. 
Also, in some cases there are very particular amplitude structures with no apparent phase correlation, such as the 
three-site structure that persist over long distances after the decomposition of the six-site in-phase configuration in 
the hexagonal geometry. It is also worthwhile to highlight the sensitive dependence of the evolution on the coupling 
parameter. For large values of the coupling, the degeneration to a single-site excitation apparently becomes more 
likely due to the fact that many of the multi-site configurations disappear through bifurcations, as has been discussed, 
e.g., in one-dimensional settings in [44J. 

There is a number of interesting directions suggested by these results for future studies. It would be relevant to 
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examine, along lines similar to the earlier work of [32[ , whether the conclusions presented herein are also qualitatively 
similar to what can be found for discrete solitons and vortices in two-dimensional Klein-Gordon chains. It would 
also be relevant to consider the continuum analogs of these results, either for cubic or for saturable nonlinearities 
(appearing, e.g., in photorefractive crystals) and see how the latter compare to the discrete theory, especially as 
concerns the qualitative structure of the linearization spectrum and the solutions predicted to be potentially stable. 
Another direction would be to explore in more detail the nonlinear dynamics, to determine the exact mechanism 
which supports the complex breathing dynamical structures, identify exact breather solutions, and to explain the 
sensitivity of the evolution on the coupling parameter. Finally, another interesting extension would be to consider 
such non-square lattices in a fully three-dimensional setting, including hexagonal-close-packed configurations. Studies 
along some of these directions are presently underway and will be reported in future publications. 

Acknowledgments 

PGK gratefully acknowledges support from NSF-DMS-0505663, NSF-DMS-0806762 and NSF-CAREER, as well 
as from the Alexander von Humboldt Foundation. The work of IK was supported by a UK EPSRC Science and 
Innovation Award. 



[1] H.S. Eisenberg, Y. Silberberg, R. Morandotti, A. R. Boyd and J. S. Aitchison, Phys. Rev. Lett. 81, 3383 (1998). 

[2] R. Morandotti, U. Peschel, J. S. Aitchison, H. S. Eisenberg and Y. Silberberg, Phys. Rev. Lett. 83, 2726-2729 (1999); H. 

S. Eisenberg, Y. Silberberg, R. Morandotti and J. S. Aitchison, Phys. Rev. Lett. 85, 1863 (2000). 
[3] D. Mandelik, R. Morandotti, J. S. Aitchison, and Y. Silberberg, Phys. Rev. Lett. 92, 093904 (2004). 

[4] D. N. Christodoulides, F. Lederer, and Y. Silberberg, Nature 424, 817-823 (2003); A. A. Sukhorukov, Yu. S. Kivshar, H. 

S. Eisenberg, and Y. Silberberg, IEEE J. Quant. Elect. 39, 31 (2003). 
[5] S. Aubry, Physica 103D, 201 (1997); S. Flach and C. R. Willis, Phys. Rep. 295, 181 (1998); P. G. Kevrekidis, K.0. 

Rasmussen and A. R. Bishop, Int. J. Mod. Phys. B 15, 2833 (2001). 
[6] N. K. Efremidis, S. Sears, D. N. Christodoulides, J. W. Fleischer, and M. Segev Phys. Rev. E 66, 46602 (2002). 
[7] J. W. Fleischer, M. Segev, N. K. Efremidis, and D. N. Christodoulides, Nature 422, 147 (2003). 

[8] J. W. Fleischer, T. Carmon, M. Segev, N. K. Efremidis, and D. N. Christodoulides, Phys. Rev. Lett. 90, 023902 (2003). 
[9] J. Yang, I. Makasyuk, A. Bezryadina, and Z. Chen, Opt. Lett. 29, 1662 (2004). 
[10] J. Yang, I. Makasyuk, A. Bezryadina, and Z. Chen, Stud. Appl. Math. 113, 389 (2004). 

[11] J. Yang, I. Makasyuk, P. G. Kevrekidis, H. Martin, B. A. Malomed, D. J. Frantzeskakis, and Z. Chen, Phys. Rev. Lett. 

94, 113902 (2005). 
[12] F. Fedele, J. Yang, and Z. Chen, Opt. Lett. 30, 1506 (2005). 

[13] D. N. Neshev, T. J. Alexander, E. A. Ostrovskaya, Yu. S. Kivshar, H. Martin, I. Makasyuk, and Z. Chen, Phys. Rev. Lett. 
92, 123903 (2004). 

[14] J. W. Fleischer, G. Bartal, O. Cohen, O. Manela, M. Segev, J. Hudock, and D. N. Christodoulides, Phys. Rev. Lett. 92, 
123904 (2004). 

[15] Y.V. Kartashov, V.A. Vysloukh and L. Torner, Phys. Rev. Lett. 93, 093904 (2004); X. Wang, Z. Chen, and P. G. Kevrekidis, 

Phys. Rev. Lett. 96, 083904 (2006). 
[16] D. Trager, R. Fischer, D.N. Neshev, A. A. Sukhorukov, C. Denz, W. Krolikowski and Yu.S. Kivshar, Optics Express 14, 

1913 (2006). 

[17] G. Bartal, O. Manela, O. Cohen, J.W. Fleischer and M. Segev, Phys. Rev. Lett. 95, 053904 (2005). 

[18] H. Trompeter, W. Krolikowski, D.N. Neshev, A.S. Desyatnikov, A. A. Sukhorukov, Yu.S. Kivshar, T. Pertsch, U. Peschel 

and F. Lederer, Phys. Rev. Lett. 96, 053903 (2006). 
[19] O. Peleg, G. Bartal. B. Freedman, O. Manela, M. Segev and D.N. Christodoulides, Phys. Rev. Lett. 98, 103901 (2007). 
[20] CR. Rosberg, D.N. Neshev, A. A. Sukhorukov, W. Krolikowski and Yu.S. Kivshar, Opt. Lett. 32, 397 (2007). 
[21] B. Freedman, G. Bartal, M. Segev, R. Lifshitz, D.N. Christodoulides and J.W. Fleischer, Nature 440, 1166 (2006). 
[22] T. Schwartz, G. Bartal, S. Fishman and M. Segev, Nature 446, 52 (2007). 

[23] J.W. Fleischer, G. Bartal, O. Cohen, T. Schwartz, O. Manela, B. Freedman, M. Segev, H. Buljan and N. K. Efremidis, 

Opt. Express 13, 1780 (2005). 
[24] Z. Chen, H. Martin, E. Eugenieva, J. Xu, and J. Yang, Opt. Express 13, 1816 (2005). 
[25] V. A. Brazhnyi and V. V. Konotop, Mod. Phys. Lett. B 18, 627 (2004). 
[26] O. Morsch and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006). 

[27] P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-Gonzalez (eds). Emergent Nonlinear Phenomena in Bose-Einstein 

Condensates: Theory and Experiment (Springer- Verlag, Heidelberg, 2008). 
[28] A. Szameit, Y. V. Kartashov, F. Dreisow, M. Heinrich, V. A. Vysloukh, T. Pertsch, S. Nolte, A. Tunnermann, F. Lederer 

and L. Torner, Opt. Lett. 33, 663 (2008). 
[29] T.J. Alexander, A.S. Desyatnikov and Yu.S. Kivshar, Opt. Lett. 32, 1293 (2007). 



10 



[30] H. Sakaguchi and B. A. Malomed. Phys. Rev. E 74, 026601 (2006). 

[31] P.G. Kevrekidis, B.A. Malomed and Yu.B. Gaididei, Phys. Rev. E 66, 016609 (2002). 

[32] V. Koukouloyannis and R.S. MacKay, J. Phys. A: Math. Gen. 38, 1021 (2005). 

[33] P K Shukla and A. A Mamun, Introduction to Dusty Plasma Physics, Institute of Physics Publishing (Bristol, 2002); 
S.V. Vladimirov, K. Ostrikov, and A. A. Samarian, Physics and Applications of Complex Plasmas, Imperial College Press 
(London, 2005). 

[34] G. E. Morfill, H. M. Thomas & M. Zuzic, Plasma crystals - a review, in Advances in Dusty Plasma Physics, Eds. P.K. 

Shukla, D.A. Mendis, & T. Desai, World Scientific (Singapore, 1997); G .E. Morfill, B.M. Annaratone, P. Bryant, A.V. 

Ivlev, H.M. Thomas, M. Zuzic and V.E. Fortov, Plasma Phys. Control. Fusion 44, B263 (2002). 
[35] O. Ishihara, announcement in 5th Int. Conf. Phys. Dusty Plasmas, Ponta Delgada, Azores (2008); private communication. 
[36] T. Misawa, N. Ohno, K. Asano, M. Sawai, S. Takamura, & P. K. Kaw, Phys. Rev. Lett. 86, 1219 (2001). 
[37] V. Koukouloyannis and I. Kourakis, Phys. Rev. E 76, 016402 (2007). 

[38] V. Koukouloyannis and I. Kourakis, Existence and stability of discrete multibreathers in hexagonal dusty plasma crystals, 
in preparation. 

[39] D.E. Pelinovsky, P.G. Kevrekidis and D.J. Frantzeskakis, Physica D 212, 20 (2005). 

[40] P.G. Kevrekidis, B.A. Malomed, Z. Chen and D.J. Frantzeskakis, Phys. Rev. E 70, 056612 (2004). 

[41] see e.g., P.G. Kevrekidis, B.A. Malomed and A.R. Bishop, J. Phys. A 34, 9615 (2001) and T. Kapitula, P.G. Kevrekidis, 

and B.A. Malomed, Phys. Rev. E 63, 036604 (2001). 
[42] D.E. Pelinovsky, P.G. Kevrekidis and D.J. Frantzeskakis, Physica D 212, 1 (2005). 
[43] T.J. Alexander, A. A. Sukhorukov, and Yu.S. Kivshar, Phys. Rev. Lett. 93, 063901 (2004). 
[44] G.L. Alfimov, V.V. Brazhnyi and V.A. Brazhnyi and V.V. Konotop, Physica D 194, 127 (2004). 



11 




FIG. 4: (Color online) The three-site configurations in the hexagonal geometry. The top two rows present the same panels as 
Fig. [2] except for the unstable unconventional case where Oi = n, and 6i = 83 = and the third is the unstable three-site 6 = 
case shown in the same format. The four panels below that display the stable (for e < 0.1) three site singly-charged vortex 
with 9 = 2n/3. They are (clockwise from top left) the modulus (for e = 0.2), phase, continuation of its principal eigenvalues as 
a function of the inter-site coupling strength e, and linear stability spectrum (for e = 0.2). Two rows are given for the 0,7r,0 
case because one of the null pairs from the AC limit becomes real (third column of the first row; the solution and its stability 
in the first and second column are shown for e = 0.025), while the other becomes imaginary (third column of the second row; 
here the solution and its stability in the first two columns are for e = 0.095). 
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FIG. 5: (Color online) The same panels as Fig. [2] except for the honeycomb geometry. The particular solutions are for 
e = 0.025 (top) and e = 0.095 (bottom). 




FIG. 6: (Color online) The same as Fig. [3]except for the honeycomb lattice. The particular solutions given are for the coupling 
constants e = 0.025 (left) and e = 0.135 (right). 




FIG. 7: (Color online) The same as Fig. [4] except for the honeycomb lattice geometry. The particular solutions shown are for 
e = 0.025 (top and third rows), e = 0.27 (second row) [the unconditionally unstable solutions], and e = 0.6 for charge 1 vortex 
solution, which is stable for e < 0.5, in the bottom set. 
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FIG. 8: (Color online) The top two rows are snapshots in the evolution of the solution corresponding to e = 0.1 from the family 
with A9 = given in the top row of Figure under a small perturbation by a random noise field to seed the instability. The 
third row shows the squared amplitude as a function of propagation distance of the relevant sites, and the inset shows a closeup 
of the small distance dynamics. Notice the structure of the robust three-site periodic structure which emerges after the original 
configuration dissolves. Below the third row panel there are two sets of images for a much larger perturbation of 25% of the 
intitial amplitude (left), where the third populated site eventually decays as well and only two sites persist for long distances, 
and a much larger coupling e = 0.3 (right), where the configuration decays very rapidly and a single site solitary wave remains. 
There is no clear correlation between the phases of either of the solutions with multiple remaining sites. 
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FIG. 9: (Color online) The same set of figures as in Fig. [8] (top row, left panels) except for the solution from the out-of-phase 
family in the bottom row of Fig. [5] are shown, also for e = 0.1. On the right a closeup of the amplitude oscillations (top) and 
phase correlation (bottom) present in the remaining breather is given. As one can see, the distance until the initial configuration 
breaks down is much longer than for the in-phase case, as expected from the linear stability analysis (cf. the inset here and in 
Fig. [8} • The ultimate surviving configuration here contains a two-site breathing structure (see bottom left panel of the squared 
amplitudes' evolution), in which there is a difference in amplitude and the phases of the the two sites are the same when the 
amplitudes are closer and opposite when they are further apart. 
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FIG. 10: (Color online) The same set as the previous images, this time for a solution from the charge-two family from the 
bottom panels of Fig. [3] with e = 0.125, i.e., large enough that a quartet of eigenvalues emerges. The long distance until 
initial breakup confirms the linear stability analysis, but a two-site (including the initially unpopulated center site) breathing 
structure persists after the disintegration of the initial structure. The inset panels show closeups of the amplitudes for shorter 
and longer propagation distances, respectively. In the right panels, one can observe that these two sites remain usually out of 
phase as their amplitudes oscillate. 



16 



2=0.1 z=3 z=6 2=10 2=20 2=100 z=202 z=300 2=402 2=500 




nnnnn nnnnn 



3 




FIG. 11: (Color online) The in-phase, three-site configuration from the third row of Fig. [4] is shown, in which a two-site 
breather persists for long propagation distances. The phases are correlated like the other unequal amplitude two-site breather 
which resulted in the evolution of the out-of-phase hexapole shown in Fig. [9] in which they become in-phase and out-of-phase 
depending on whether their amplitudes are similar or considerably different, respectively. 
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FIG. 12: (Color online) The evolution of the in-phase six-site configuration with the honeycomb lattice geometry from the 
top row of Fig. [5] is given above (left) for e — 0.1 and (right) for e = 0.3. As in the hexagonal case shown in Fig. [8] for 
e = 0.1 a multi site structure persists over a long distance, although now it is comprised of four sites, two pairs of out-of-phase 
breathers with comparable amplitude (the phase structure is not shown, but each pair is comparable with that of Fig. I10|l . 
This interesting difference inspired us to continue the dynamical evolution for a longer distance, and the structure did indeed 
persist up to another order of magnitude. Even with the much larger perturbation of 25% of the initial amplitude (not shown) 
as opposed to 5%, a very similar four site structure persists for a long distance, although the degeneration of the other two 
sites is very rapid. The same robustness to perturbation is found for e = 0.3, although a two site unequal amplitude breather 
remains, which oscillates between in-phase and out-of-phase (not shown, but same as the unequal amplitude breather in Fig. 
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FIG. 13: (Color online) The six-site doubly-charged honeycomb lattice vortex for e = 0.135 from the bottom set of panels in 
Figure[6]is significantly more stable than the singly-charged counterpart. All original sites remain populated for a long distance, 
up to z = 400, and, when the two sites eventually decrease in amplitude, the remaining four reshape into a four-site breather. 
Two of the sites remain close in amplitude and out-of-phase, while one has larger and the other has smaller amplitude and 
these oscillate between in-phase and out-of-phase in the same manner as the others, such as those in Fig. [9] (see panels on 
the right). The inset features a closeup image of the complex oscillations of the four sites. At a very long distance, close to 
z = 2000, they reshape in amplitude and phase, becoming two pairs of out-of-phase breathers, like the in-phase hexapole from 
Fig. \T%\ ultimately does, although the dynamics is not followed further to see if this structure persists. 
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FIG. 14: (Color online) The dynamics for the 0, ir, honeycomb lattice configuration from the top rows of Fig. [7] for s = 0.1 
This solution persists for very long propagation distances despite the linear instability. Moreover, the relative phase structure 
persists, although the one site that is out-of-phase with the other two oscillates from one to another among the three (see the 
right panels). 



